Charge constrained density functional molecular dynamics for simulation of 
condensed phase electron transfer reactions 
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We present a plane-wave basis set implementation of charge constrained density functional molec- 
ular dynamics (CDFT-MD) for simulation of electron transfer reactions in condensed phase systems. 
Following earlier work of Wu et al. Phys. Rev. A 72, 024502 (2005), the density functional is mini- 
mized under the constraint that the charge difference between donor and acceptor is equal to a given 
value. The classical ion dynamics is propagated on the Born-Oppenheimer surface of the charge 
constrained state. We investigate the dependence of the constrained energy and of the energy gap 
on the definition of the charge, and present expressions for the constraint forces. The method is 
applied to the Ru'^^-Ru^^ electron self-exchange reaction in aqueous solution. Sampling the vertical 
energy gap along CDFT-MD trajectories, and correcting for finite size efi'ects, a reorganization free 
energy of 1.6 eV is obtained. This is 0.1-0.2 eV lower than a previous estimate based on a contirmum 
model for solvation. The smaller value for reorganization free energy can be explained by the fact 
that the Ru-0 distances of the divalent and trivalent Ru-hexahydrates are predicted to be more 
similar in the electron transfer complex than for the separated aqua-ions. 



I. INTRODUCTION 

Atomistic simulation of electron transfer reactions has 
become a major task in computational physics and 
chemistry. This field was pioneered in the Eighties 
by the works of Warshel0, % 0], Kuharski et a/[i, Q 
and others [1, 0, [I]- Electron donor and acceptor and 
the solvent were modelled with classical potential en- 
ergy functions and the diabatic free energy profiles were 
sampled with (biased) molecular dynamics simulation. 
These early simulations gave unprecedented insight into 
electron transfer reactions and provided a first numer- 
ical confirmation for the validity of the linear response 
approximation of the Fe^+-Fe'^~'" electron self-exchange 
reaction thus confirming a crucial assumption in Mar- 
cus theory of electron transfer 0. 

The application of modern density functional molecu- 
lar dynamics to electron transfer reactions has not been 
successful until recently [13] ■ The reason is that common 
exchange-correlation functionals are of limited use for 
this task due to their tendency to erroneously delocal- 
ize electrons, thus prohibiting an accurate modelling of 
charge transfer states. This deficiency termed electron 
self-interaction or delocalization error [IT, is intrin- 
sic to GGA and hybrid density functionals unless special 
care is taken in their parametrization. fn parallel to the 
development of functionals with minimal self-interaction 
error jl^ J^], a number of correction schemes have been 
proposed to minimize the delocalization error of exis ting 
and coinputationally inexpensive density functionals 1 1 5l. 
m, [13, [3 [H, including DFT-|-U;2a 2l[, [H, formula- 
tion of a penalty density functiona llld l , a nd constrained 
DFT[21, [13, [H, [H, [13, M [H [13, [3l[[33, [H . Develop- 
ment of such schemes are particularly important in den- 



* Electronic address: 'j b376@cam.ac.uk] 



sity functional molecular dynamics, where efficient com- 
putation of the exchange-correlation energy and forces is 
absolutely crucial. 

In recent years several groupsflO, [H, [11] have greatly 
advanced the development of the constrained DFT ap- 
proach. In this method an external potential is added to 
the Kohn-Sham equations preventing an excess electron 
or hole from wrongly delocalizing over donor and accep- 
tor ions. The external potential is varied until a given 
constraint on the density, for instance a charge or spin 
constraint, is satisfied. The search for this external po- 
tential, usually carried out within a Lagrange multiplier 
scheme, introduces a second iteration loop in addition 
to the usual self-consistent iteration of the Kohn-Sham 
equations. Due to the localized (or diabatic) nature of 
the constrained states the delocalization error of common 
density functionals is reduced, though not eliminated. 
Hence one can expect that common density functionals 
describe constrained states as well as any other states 
where delocalizing electrons are not present. 

The construction of charge localized states seems to be 
a somewhat artificial procedure since the external poten- 
tial and charge are all but uniquely defined. For large 
donor-acceptor distances one can expect, however, that 
the details of charge definition are less relevant. Thus 
it is appealing to use CDFT to describe the charge lo- 
calized (diabatic) states of long-range electron transfer 
reactions. However, as it is possible to extract approxi- 
mate electronic coupling matrix elements between CDFT 
states ^26,], one can construct approximate Hamiltonian 
matrices in the space of constrained states (CDFT con- 
figuration interaction) . The adiabatic states diago- 
nalizing this Hamiltonian suffer less from the delocaliza- 
tion error than unconstrained states and can be used to 
describe short ran ge p henomena such as chemical bond 
break reactions p9l,T30| . 

All of the above mentioned calculations (except the 
ones described in Refs. [HI, [1^) were carried out in 



2 



the gas-phase using quantum chemistry codes. In this 
work we present a plane-wave basis set implementation 
of charge constrained density functional molecular dy- 
namics (CDFT-MD) in the Car-Parrinello molecular dy- 
namics code (CPMD)(3^]. This allows us to study the 
electron self-exchange of the Ru^+-Ru''+ ion pair in aque- 
ous solution at finite temperature with the solute and 
solvent treated at the same density functional level of 
theory. To our best knowledge only one such simulation 
has been reported previously by Sit et al. , who devel- 
oped a penalty density functional approach to simulate 
the isoelectronic Fe^+-Fe^+ electron self-exchange reac- 
tion in aqueous solution. A previous attempt to use a 
charge restraint rather than a charge constraint to drive 
a charge transfer reactions with Car-Parrinello molecular 
dynamics was also reported [35]. However, it turned out 
that at reasonably high restraining forces the charge re- 
straint does not provide enough bias for obtaining charge 
localized states. 



In the present work we define the constraint as 
the charge difference between the electron donating 
Ru^+-hexahydrate and the electron accepting Ru'^+- 
hexahydrate. The charge constrained state is energy op- 
timized at each molecular dynamics step and the dynam- 
ics of the system propagated on the constrained Born- 
Oppenheimer surface. We then construct Marcus-type 
free energy profiles by calculation of the vertical energy 
gap between the two charge transfer states along the 
molecular dynamics trajectory. The reorganization free 
energy obtained, about 1.6 eV after correction for finite 
size effects, is in fair agreement with experimental esti- 
mates and other computational studies indicating that 
CDFT-MD can be successfully applied to electron trans- 
fer problems in the condensed phase. 



This paper is organized as follows: First we briefly re- 
view the theoretical background of CDFT. Then we ad- 
dress the dependency of CDFT energies on the choice of 
the constraining potential for gas phase and condensed 
phase electron transfer systems. Thereafter our molec- 
ular dynamics implementation of CDFT is validated by 
finite temperature simulation of the molecule. The 
main objective of the present work, the charge con- 
strained density functional MD simulation of the aqueous 
Ru'^+-Ru'^+ electron self exchange reaction is presented 
thereafter. The results are discussed in light of experi- 
mental results, and previous classical molecular dynamics 
and continuum studies. At the end of this work we give 
a perspective on future work planned. In the appendix 
we give explicit expressions for constraint forces due to 
the charge constraint and summarize relevant technical 
details of CDFT-MD calculations. 



II. THEORY 

A. Charge constrained density functional 
molecular dynamics 

CDFT and its working equations have been presented 
previously in a number of papers [H, [l^, [H, H^, 
Isil . [3^ . [331 . Here we follow the work of van Voorhis and 
co-workers [2^ and give a short summary of the equations 
pertinent to our molecular dynamics implementation. 

In CDFT the usual energy functional E[p] is minimized 
with respect to the electron density p under the condition 
that the sclcronomic constraint 



w{v)p{v) dv = Nc, 



(1) 



is satisfied, where Nc is a real number. The weight 
function w{r) on the left hand side of Eq. [T] defines the 
constraint, for instance the charge of an atom, molecule 
or molecular fragment, or the charge difference between 
groups of atoms. Nc defines the value of the constraint. 
Both quantities are input parameter that remain un- 
changed during the constrained minimization. Using the 
Lagrange multiplier technique the new energy functional 
to be minimized is given by 

W[p, V] = E[p] + V( j w{y)p{v) dr ~ N^), (2) 

where V is an as yet undetermined Lagrange multiplier. 
In order to find the minimum of W[p, V] with respect to 
(wrt) p and V, W is minimized with respect to p for a 
given V, the gradient 



SW 
W 



w{r)p{r) dr — Nc 



(3) 



at the minimizing density is determined to generate the 
next iteration step for V, W is minimized for the new 
value of V wrt p and this procedure is repeated self- 
consistently until convergence is reached, that is when 
SW/SV\p = and SW/6p\v = 0. Eflicient Newton meth- 
ods can be employed for this minimization procedure 
since second derivatives of W wrt V are available an- 
alytically via density functional perturbation theory [25j 
or numerically as finite differences of Eq. [3l For practical 
applications we define in addition to the usual wavefunc- 
tion convergence criterion a second convergence parame- 
ter for the constraint. 



C > 



w{r)p{r) dr — Nc 



(4) 



C is a measure of how accurately the charge constraint 
Eq.[T]is fulfilled. 

It is interesting to note that Eq. ([2]) implies the fol- 
lowing exact relation for the energy difference between 
two constrained states A and B, AEab, and the continu- 
ous set of Lagrange multipliers connecting the two states. 
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centered weight function Wi 



AEab 



(5) 



The Lagrange multipher can thus be interpreted as the 
force along the charge coordinate iVc. The paraUel to 
thermodynamic integration is intriguing. 

CDFT-MD simulation or constrained geometry opti- 
mization can be carried out provided one takes into ac- 
count the additional forces Fd arising from the constraint 
term on the right hand side of Eq. [21 Adopting the 
Hellmann-Feynman theorem the total force on atom i 
at position is given by 



where Ft, 



F,, 



'dW/dRi, F, = -dE/d'R^ and 



dw{r,R) 



(6) 



(7) 



Note that the weight function can depend on the coordi- 
nates R of all atoms in the system. An explicit expres- 
sion for the derivative in the integrand of Eq. [7] using the 
weight function defined below is given in appendix |B] In 
molecular dynamics simulation the Lagrange multiplier 
V, the total energy and forces have to be calculated at 
each time step. The computational bottleneck is the it- 
erative search for V. To improve the efficiency we have 
implemented an extrapolation scheme for the Lagrange 
multiplier using Lagrange polynomials. For details we 
refer to appendix ICl 



B. Definition of the charge constraint 

The charge constraint is fully defined by the weight 
function w in the integrand of Eq. [T] and the actual value 
of the constraint, Nc. In principle there are an infi- 
nite number of ways of how to choose the weight func- 
tion. In practice one chooses the weight so that Eq. [1] 
corresponds to some common charge definition. One is 
then left to investigate how much the results depend on 
the charge definition used. In previous work MuUiken, 
Loewdin and a Becke real space integration scheme have 
been testedfl^. While for short donor-acceptor distances 
the energy of the constrained state was strongly depen- 
dent on the weight used, for medium to large distances 
the dependence on the weight was reasonably small. The 
real space density integration scheme was found to give 
best overall performance ^26, .29] . 

Building on this previous work we use a slightly dif- 
ferent real space density integration scheme for charges, 
the one according to Hirshfeldjs^. The Hirshfeld charge 
Qi of an atom i at position R^ is obtained by integration 
of the total electron density p multiplied with an atom 



Qi = Zi - J Wi{r,R)p{r) dr 



Wi(r, R) = 



Pi{r - Rj) 
EliPj(r-Rj) 



(8) 
(9) 



where Zi is the core charge in pseudopotential calcula- 
tions (or the charge of the nucleus in all-electron calcula- 
tions), and N is the total number of atoms. The weight 
function is constructed from the unperturbed promolec- 
ular densities of atoms i, pi, 



(r-R,) = A(r) = ^n,|V4(r)P 



(10) 



where r= |r — R^j. The sum ranges over the radial part 
of the promolecular atomic orbitals ipl {r) with occupa- 
tion number Uj. The weight Wi is close to unity up to a 
distance of about one atomic radius and goes to zero ac- 
cording to the decay of the promolecular atomic density. 
Note that the sum of the Hirshfeld charges of all atoms 
is equal to the total charge of the system. 

A natural choice for the constraint for electron trans- 
fer reactions is the charge difference between a set D of 
atoms comprising the electron donor and a set A of atoms 
comprising the electron acceptor. 



C- 



where 



w 



WD - WA 



Ef=iPj(r-Rj) 



(12) 

and Pi is given by Eq. [TOl The constant C in Eq. [TT] is 
equal to the difference in the core charges, C = X]ieD ^i" 
J2ieA merely causing a constant shift of W (Eq. ^ 
by VC. The weight function Eq. [12] used for the simu- 
lation of the Ru^^-Ru'^"'" electron self-exchange reaction 
in aqueous solution is illustrated in Fig. [TJ The donor 
(acceptor) atoms are comprised of Ru^"*" (Ru'^+) and all 
atoms of the first solvation shell of Ru^"*" (Ru'^'*') . Note 
that the sign of the weight function changes sharply at 
the interface of the donor-acceptor complex. Assuming 
that a charge equivalent to one electron is transferred, 
the constraint value A^c is set equal to 1 for the initial 
state and equal to -1 for the final state. 

While relying on the basic charge definition Eqs. [SHU 
we have investigated different functional forms for the 
densities pi that define the atomic weight function Wi of 
Eq. [9] The results are presented in section IIV CI Techni- 
cal details concerning the calculation of the weight func- 
tion can be found in appendix [XI 



C. Electron transfer free energy curves 

In Marcus theory 9] electron transfer reactions are de- 
scribed by two charge localized or diabatic free energy 
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FIG. 1: Isosurfaces of the weight function Eq. [T2]for values 
w = 0.6 (red) and w — —0.6 (blue). The snapshot is taken 
from a charge constrained MD simulation of the Ru^"^-Ru'^^ 
electron transfer complex in aqueous solution. Ru ions are 
depicted as green spheres and first shell water molecules are 
depicted in ball and stick representation. All atoms compos- 
ing the Ru-hexahydrate in the upper right (lower left) part 
of the figure are donor (acceptor) atoms. Solvent molecules 
not included in the charge constraint are depicted in stick 
representation. 



curves, one for the reactant state (A) and one for the 
product state (B). Electron transfer is assumed to occur 
at the crossing point of the two curves. In the limit of lin- 
ear response the two free energy curves are parabolic and 
the activation free energy for electron transfer is given by 



4A 



(13) 



where is the reaction free energy and A is the reorga- 
nization free energy. The latter is inversely proportional 
to the curvature of the free energy parabola and is a mea- 
sure for the free energy required to distort the equilibrium 
configuration of one diabatic state to the equilibrium con- 
figuration of the other diabatic state while staying on the 
same free energy curve. For electron self-exchange reac- 
tions AA = and the activation free energy is entirely 
determined by the reorganization free energy. 

Following WarshclfT] the diabatic free energy curves 
defining the reorganization free energy can be obtained 
by sampling the vertical energy gap AE 



AE(R) =^;b(R)--Ba(R) 



(14) 



using molecular dynamics simulation. Ea (R) and Eb (R) 
are the charge localized (diabatic) potential energy sur- 
faces and the vector R denotes the coordinates of all 
atoms in the system. The relative Landau free energy 
AnfiAE) along this coordinate is given for state M, 
M = A, B, by 



where fce is the Boltzmann constant, T the temperature 
and pm{AE) the probability distribution of the reac- 
tion coordinate in state M. If the free energy curves 
are parabolic the reorganization free energy is equal to 
the average vertical energy gap. 



A = (AE), 



(16) 



where (•••)a denotes the usual canonical average for 
state A. 

The free energy curves Eq. [15] can be obtained by 
sampling configurations with charge constrained density 
functional molecular dynamics in state A as described 
above, followed by calculation of the vertical energy gap 
Eg. [T4l between the constrained states A {Nc = l, Eq. fTTj) 
and B {N^ = — 1) for the set of sampled configurations. 
However, unbiased equilibrium simulations give only ac- 
curate results close to the free energy minimum of Am 
and are of limited use for regions of AE far away from the 
minimum. Fortunately, due to the exact linear free en- 
ergy relation between the free energy gap and the energy 

gapSaSii, 



Ab{AE) - AAiAE) ^ AE 



(17) 



Am{AE) ^ -kBTlnpM{AE), 



(15) 



it is possible to calculate a good part of the curve at 
high free energies accurately from equilibrium simula- 
tions. Thus, using information from two distinct regions 
of AE one can construct a reasonably accurate free en- 
ergy profile without the use of computationally expensive 
enhanced sampling methods [sR [39| . 



III. COMPUTATIONAL DETAILS 

Charge constrained density functional molecular dy- 
namics has been implemented in the CPMD code[3J|. 
Unless stated otherwise, all calculations were carried 
out with the BLYP |40s ^] functional using a reciprocal 
space cutoff of 70 Ry, TrouUier-Martins pseudopotentials 
for nuclei-|-core electrons [i^ , pseudoatomic densities for 
construction of the weight function Eq. I12[ a convergence 
criterion for the wavefunction gradient of 1 x 10~^H and a 
charge constraint convergence criterion defined in Eq. (|4]) 
of C = 5 X lO^^e. All calculations were carried out in the 
lowest spin state. 

Ru^+-Ru^+ electron self-exchange was simulated in a 
periodic box of dimension 14.5 x 11.35 x 11.35A con- 
taining two Ru ions and 63 water molecules. The donor 
and acceptor groups are comprised of the ion and the 
six water molecules forming the first solvation shell, 
Ru(H20)^+ and Ru(H20)^+, respectively, see Figured] 
The charge constraint Eq. [11] was set equal to 1 corre- 
sponding to the reactant diabatic state A. Pseudoatomic 
densities for construction of the weight function Eq. [T2l 
were used, i.e. the charges correspond to the definition 
of Hirshfeld[36j. The Ru^"*" and Ru'^"'" aqua ions are both 
low spin as opposed to the corresponding aqua-ions of 
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Fe2+ and Fe^+fT^. Hence, the dublet state was cho- 
sen for simulation of the aqueous ET complex. The ini- 
tial configuration was taken from an equilibrated classical 
molecular dynamics trajectory carried out in a previous 
investigation of the same reaction 43[ . The distance be- 
tween the two Ru ions was fixed at 5.5 A using the 
RATTLE algorithm The system was simulated in 

the NVT ensemble at 300 K using a chain of Nose-Hoover 
thermostats [45I of length 4 with a frequency 1000 cm~^. 
To increase the efficiency of charge constrained Born- 
Oppenheimer molecular dynamics we used an extrapo- 
lation scheme for prediction of the Lagrange multiplier 
V as described in appendix [Cl an extrapolation scheme 
for the initial guess of the wavefunction, slightly higher 
convergence criteria of C = 5 x 10~^ e for the constraint 
convergence and of 2 x 10~^ H for the wavefunction gra- 
dient, and an MD time step of 40 au — 0.96 fs. With 
this setup the average number of iterations for V were 
~ 2 — 3 per molecular dynamics time step. Accordingly, 
the computational overhead compared to standard Born- 
Oppenheimer molecular dynamics without charge con- 
straint is about a factor of 2-3. The average drift of 
the conserved energy was —9.7 x 10~^H/atom/ps along 
a trajectory of length 6.6 ps. This is somewhat large 
but still acceptable for our purposes. Better energy con- 
servation can be obtained if tighter convergence criteria 
and a smaller time step is used. For calculation of ther- 
mal averages the first ps of dynamics was discarded. The 
energy gap Eq. [14] was calculated for 225 equidistantly 
spaced snapshots taken from the last 5.5 ps of the trajec- 
tory using the same convergence criteria as for the MD 
simulation. For calculation of the product diabatic state 
B the constraint Eg. [TT] was set equal to -1. 



IV. RESULTS AND DISCUSSION 

Before we present the results for electron self-exchange 
in the condensed phase we report on a series of test calcu- 
lations carried out for simple electron transfer systems in 
the gas phase. These include a basic validation of charge 
constrained single point calculations and molecular dy- 
namics, calculation of the dissociation curve of HJ and 
HeJ to test the correct long-range behaviour of CDFT, 
and an investigation of the dependence of the constrained 
state energies on the weight function used. 



A. (CO)- 

As a first test of our implementation we carried out a 
series of charge constrained wavefunction optimizations 
on the (C0)~ molecule. The distance between the C and 
the O atom was chosen to be 2A in order to avoid the 
problem of spin contamination which occurs for larger 
distances. The oxygen atom is chosen as the electron 
donor and the carbon atom as the electron acceptor. The 
excess electron is transferred from the oxygen to the car- 



bon atom by increasing the charge difference A^c (Eq. fTTj) 
between the atoms in small steps from -1 to 1. The La- 
grange multiplier V and the potential energy E of the 
constrained states are shown in Fig. [5] as a function of 
N,. 




FIG. 2: Lagrange multiplier V (A) and potential energy E 
(B) versus constrained charge difference Ac between C and O 
atom in (CO)~. In (B) the computed energies are indicated 
by circles and the running integral - V dNc + E[C'0] 
by straight lines. The potential energy and charge difference 
for unconstrained wavefunction optimization of (C0)~ is in- 
dicated by a cross. 

According to Eq.[5]the energy difference between C^O 
and CO^ should be equal to the negative of the integral 
of V over A^c- Indeed, the difference between the end- 
points is AEab = E[C~0]-E[CO~]=9m8 mH whereas 
integration gives — J^^ V dNc = 9.084: mH. The numbers 
match within the given convergence criteria. Further- 
more, the potential energy obtained from unconstrained 
calculations (denoted by a cross in Fig. [2]) lies on the 
constrained potential energy curve at A'c ~ 0. Thus, the 
unconstrained state can be reproduced by constrained 
wavefunction optimization if a constraint value is used 
that is equal to the charge difference of the unconstrained 
state. 



B. Dissociation of and Hej 

The wrong dissociation curve for H^ is probably one of 
the most spectacular failures of GGA and hybrid de nsity 
functionals. The reason for this failure is well known |ll|. 
Due to the wrong scaling behaviour of these functionals 
wrt electron number the charge delocalized state is pre- 
dicted to be lower in energy than the charge localized 
state, even though the two states are degenerate in the 
limit of large inter-nuclear separation distance in exact 
theory. Using CDFT configuration interaction Wu et al. 
could circumvent this problem [2^. The authors reported 
a dissociation curve that matched the exact Hartree-Fock 
curve remarkably well at the equilibrium region and at 
long range. 
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Here we calculate the dissociation curve for and 
HeJ using a single charge constrained state in order to 
test the correct long range behaviour of our CDFT im- 
plementation. The constraint is again the charge differ- 
ence between the two atoms. In case of HeJ iVc was set 
equal to 1 for all distances. For the constraint value 
was set equal to the charge difference obtained when the 
promolecular pseudoatomic reference orbitals were used 
for construction of the initial wavefunction(29| . This was 
necessary because a value N^ — l would only be obtained 
at very high external potentials for which the wavcfunc- 
tion would not converge. Thus, the constraint value for 
changes with distance from about 0.86 at 2A to 1.0 
at 4.4A and larger distances. 

The energies of the charge constrained states of 
and HeJ relative to the (unconstrained) energy of the 
isolated fragments are shown in Fig. [31 together with the 
exact Hartree Fock curve for and the essentially exact 
FCI curve for He J ^] . The CDFT curves agree well with 




FIG. 3: Dissociation curves for (A) and He J (B). Con- 
strained and unconstrained DFT results are indicated by 
squares and circles, respectively. The exact Hartree-Fock 
curve for and the curve for He J at the FCI/aug-cc-pV5Z 
level of theory[4^ are indicated by triangles, respectively. 

the exact curves for distances larger than SA. This is a 
considerable improvement with respect to the results of 
unconstrained optimizations, which give the wrong dis- 
sociation limit. It also shows that the optimization in 
the external potential yields the correct electronic state 
of the unconstrained isolated fragments. The significant 
deviation for smaller distances is due to the fact that 
only a single charge localized state is used, which can 
not describe chemical bonding between the fragments. 



The deficiency at short range can be cured using the 
CDFT-configuration interaction approach introduced in 
Ref. fH. 



C. Dependence on the weight function 

A crucial issue in CDFT calculations is the depen- 
dence of the numerical results on the definition of the 
charge constraint. While relying on the basic real-space 
charge definition Eqs. [HEl we have investigated differ- 
ent functional forms for the densities pi that define the 
atomic weight function Wi of Eq. [S] Besides numerical 
pseudoatomic densities we have tested a minimal basis 
of Slater functions and a minimal basis of Gaussian func- 
tions, 



Pi{r - Rj) = 



(18) 



where ak is the width of particle species k and Nci,i de- 
notes the number of electrons of the isolated atom i. Us- 
ing Gaussian functions one can vary the decay of the 
weight functions systematically simply by changing the 
exponent of the densities in Eq. [THl For slowly decay- 
ing densities the weight function changes sign smoothly, 
while for densities that match the pseudoatomic densities 
the sign changes sharply. In the limit of an infinitely fast 
decaying density pi the electron density p at a given point 
in space is assigned to the atom closest to this point. 
This corresponds to the charge definition according to 
Voronoi. For the real space integration of Eq. [8]we have 
introduced a radial cutoff Rc (see appendix lA)) . Its influ- 
ence on the constrained energy is also reported here. The 
charge constrained energies of ZnJ obtained for different 
weight functions are summarized in table HI The Zn-Zn 
distance was chosen to be r = 4A and r — 2.5 A. The 
latter is equal to twice the covalent radius of Zn. The 
constraint was again the charge difference, A'c = l. Con- 
sidering hole transfer at a distance r = 4A, we find that 
the change in energy is only a few mH when the width 
of the Gaussian function is varied from 0.5 to 1.0 A. If 
larger values for the widths are used the weight function 
becomes unphysically smooth and the energy increases. 
The energy obtained with pseudoatomic functions dif- 
fer by not more than 2 mH from the energies obtained 
with Gaussian functions. A somewhat larger deviation 
is obtained for Slater functions. The cutoff i?c = 3.762A 
for truncation of the weight function is sufficient for all 
weight functions. Overall, the dependence of the results 
on the weight function used is reasonably small. This 
is not the case for hole transfer at a smaller distance of 
r — 2.5 A. Here the details of the weight function are im- 
portant. The variation in energy are a few ten mH, an 
order of magnitude larger than at r = 4A. Hence, our 
results indicate that charge constrained states are well 
defined only if the distance between donor and acceptor 
is larger than at least the sum of their covalent radii. 
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TABLE I: Charge constrained energies for hole transfer in 
ZnJ. Pseudoatomic, Slater and Gaussian denote the func- 
tional form of the densities pi (Eq. I10|l that define the atomic 
weight function Wi of Eq. [S] cr is the width of the Gaussian 
function Eq. [18] and i?c denotes the radial weight cutoff (see 
appendix O . 



r = 4.0A 


WPl 0'Vlf 




R fAl 


j^iiuigy [iij 


pseudoatomic 




3.762 
5.121 


-113.7834 
-113.7834 


Slater 




3.762 
5.121 


-113.7780 
-113.7783 




0.5 


3.762 


-113.7859 




0.732 


3.762 


-113.7849 


Gaussian 


1.0 
2.0 


3.762 
3.762 


-113.7818 
-113.7633 




1.0 


2.56 


-113.7813 




1.0 


3.0 


-113.7814 




1.0 


5.0 


-113.7819 




r = 2.5A 


weight 


a [A] 


Rc [A] 


Energy [H] 


pseudoatomic 




3.762 
5.121 


-113.7744 
-113.7745 


Slater 




3.762 
5.121 


-113.7471 
-113.7473 




0.5 


3.762 


-113.7824 




0.732 


3.762 


-113.7755 


Gaussian 


1.0 
2.0 


3.762 
3.762 


-113.7613 
-113.6495 




1.0 


2.56 


-113.7781 




1.0 


3.0 


-113.7612 




1.0 


5.0 


-113.7617 



As we are primarily interested in the aqueous Ru'^+- 
Ru^+ electron self-exchange reaction (see section IIV E[) 
we have investigated the dependence of the vertical en- 
ergy gap Eq. [13] of this system on the functional form 
of the weight used. For this purpose we have taken a 
snapshot from the constrained MD simulation of aqueous 
Ru^+-Ru'^"'' and calculated the two charge constrained 
states Ea (7Vc = -1) and Eb (iVc = -l). For the Gaus- 
sian weight function cth = cq = 0.6 A and ctru = 1-0 A is 
used. See section Hill for further details. The results are 
summarized in table |TT1 The energy gaps differ by less 
than 0.18 eV for the three weights considered and are 
within 0.04 eV for pseudoatomic and Gaussian weight 
functions. This variation is not insignificant and should 
be considered as a lower limit of the error of the results 
presented in section IIV El More than the gap energies 
varies the charge of the electron donor Ru^+(H2 0)6 (re- 



TABLE II: Dependence of the energy gap Eq.[T4]on the weight 
function for the aqueous Ru^^-Ru''''' electron self-exchange 
reaction. The charge is the sum of the atomic charges of 
Ru=^+(H2 0)6 in the final ET state. 



weight 


energy gap [eV] 


charge [e] 


pseudoatomic 


1.587 


1.5663 


Slater 


1.732 


1.2639 


Gaussian 


1.552 


1.7764 



call that only the charge difference between donor and 
acceptor is constrained). The sensitivity of the results 
on the weight used is probably a consequence of the close 
approach of the first shell water molecules that bridge the 
two Ru-ions. These water molecules form strong hydro- 
gen bonds that make the constrained energy and charges 
susceptible to details of the weight function at the inter- 
face of the two Ru-complexes. 



D. CDFT-MD for 

The molecular dynamics implementation of CDFT is 
tested by simulating an isolated molecule on the 
Born-Oppenheimer surface of a single charge-constrained 
state. A charge difference Nc — 0.5 e is enforced giving 
an average charge of 0.25 e respectively 0.75 e on the two 
H atoms. The system is simulated in the NVE ensem- 
ble at a temperature of approximately 300 K using the 
Velocity Verlet algorithm. A value of 5 x 10^® H for the 
convergence of the wavefunction gradient is used and a 
timestep of 20au ~ 0.48 fs. In order to assess the de- 
pendence of energy conservation on the convergence cri- 
terion for the charge constraint, Eq. 21 we calculated a 
series of trajectories of length 1 ps for different values 
of C. The total linear drift of the conserved energy as a 
function of C is shown in Fig. 01 The drift is less than 

0.04 n ' ' 




. . . , , 

le-06 le-05 0.0001 



Constraint convergence [e] 

FIG. 4: CDFT-MD simulation of H+. The linear drift in 
the conserved energy is shown as a function of the charge 
constraint convergence criterion C defined in Eq. [J] 
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1 X 10~^ mH/atom/ps at C = 1 x 10~^e showing that 
the total energy is essentially conserved if a tight con- 
vergence criterion for the charge constraint is applied. 
The sharp rise of the total energy drift at a value of 
C = 1 X 10~^ e is due to the fact that the system essen- 
tially behaves like a harmonic oscillator. Small deviations 
from the target constraint value can lead to a resonance 
between the bond vibration and the external constraint 
potential. This can cause instabilities in the integration 
of the equations of motion, which can even lead to dis- 
sociation of the molecule. Fortunately, this behaviour is 
rather exceptional. The resonance effect is dampened in 
larger systems allowing us to use less strict convergence 
criteria than for H J . 



E. Ru "''-Ru + electron self-exchange in aqueous 
solution 

We finally present our results for the CDFT-MD sim- 
ulation of Ru^+-Ru'^"'" electron self-exchange in the con- 
densed, aqueous phase. The charge constraint is chosen 
as the charge difference between the electron donating 
group, Ru^+(H20)6, and the electron accepting group, 
Ru'^+(H2 0)6, and the constraint value is A^c = 1- This 
choice is motivated by the fact that one wants to trans- 
fer a charge corresponding to one electron from the donor 
to the acceptor complex. The six water molecules form- 
ing the first coordination shell are included in the con- 
straint as the redox active orbitals are delocalized over 
the metal and the first coordination shell. The distance 
of the two Ru ions is constrained to 5.5 A. The same dis- 
tance was used in a previous classical molecular dynamics 
simulation of this reactionji^. No other constraints on 
the dynamics are imposed. Details for the CDFT-MD 
simulation are summarized in section [TTTl 

During CDFT-MD the highest occupied majority spin 
orbital (HOMO) of the FT complex is correctly located 
on the donor complex, and the lowest unoccupied mi- 
nority spin orbital (LUMO) is located on the acceptor 
complex. As expected, the two molecular orbitals are 
composed of a d orbital of the metal t2g manifold and the 
p orbitals of the ligands. Similarly to the separated aqua 
ions, there is no mixing with orbitals of solvent molecules 
beyond the first solvation shell. 

As the charge difference is constrained, only, the abso- 
lute charge of donor and acceptor complex is free to vary 
during the dynamics run. The charge fluctuations are 
very small, however, a = 0.05 e. The average charge of 
the Ru2+(H20)6 complex, 0.52e, and of the Ru3+(H20)6 
complex, 1.52 e, are significantly smaller than their for- 
mal charges of -1-2 e and -1-3 e, respectively. They are, 
however, similar to the charge of a single Ru^+(H20)6 
(Ru^+(H20)6) ion in aqueous solution, 0.75 e (1.15 e). 
Thus, the charge constraint localizes an excess charge of 
— 0.23e (0.37 e) on the donor (acceptor) complex relative 
to the charge of the isolated aqua ions. The charge of the 
remaining solvent, 2.96 e, is more than half of the total 



system charge. Although some charge transfer between 
the ET complex and the solvent is expected, the mag- 
nitude of this effect seems rather large. The reason for 
this is not clear, but one may speculate that the BLYP 
functional tends to delocalize the total system charge, an 
effect that might be enhanced when periodic boundary 
conditions are used for simulation of the aqueous phase. 
Yet, the large magnitude of charge transfer to the solvent 
is not a particular feature of CDFT, because this effect 
already occurs for standard (unconstrained) GGA-DFT 
calculations on solutions containing a single ion. 

In order to assess the effects of the charge constraint 
on the coordination geometry, we calculated the metal- 
oxygen radial distribution functions (5ruo('')) of Ru^+ 
and Ru^+ in the electron transfer complex and compare 
with the radial distribution function of the single aque- 
ous ions as obtained from standard (unconstrained) Car- 
Parrinello molecular dynamics simulation [47l. The 
result is illustrated in Fig. [5l In the unconstrained sim- 




r[A] 



FIG. 5: First peak of the radial distribution functions of Ru^^ 
and Ru"^''' with the oxygen atoms of the solvent molecules. 
Curves labeled 'CDFT' were obtained from present CDFT- 
MD simulations of the solvated ET complex, curves labeled 
'single' were obtained from standard Car-Parrinello MD sim- 
ulation of a single ion in aqueous solution 47, 48], and curves 
labeled 'unconstrained' were obtained from standard Car- 
Parrinello MD simulation of the solvated ET complex. Note 
that the curves for Ru^+ and Ru'^"'' coincide for the latter sim- 
ulations. All distribution functions were smoothed by convo- 
lution with a Gaussian of width O.OSA. 

ulations of the single aqua ions the average Ru-0 bond 
distances are 2.17 A for Ru^+ and 2.09 A for Ru^+jiS,!!!. 
The difference in distance is significantly smaller in the 
ET complex, 2.15 A for Ru'^+ and 2. 13 A for Ru^+ . The 
deformation of the two complexes to a more similar co- 
ordination geometry must be attributed their strong in- 
teractions at a rather short Ru-Ru distance of 5.5 A. The 
solvation shells of the two ions interpenetrate and the wa- 
ter molecules bridging the two Ru ions form 1 — 2 strong 
hydrogen bonds during the course of the simulation. As 
expected, the solvation structure of the two ions in the 
ET complex is virtually identical if the charge constraint 
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is not imposed. The radial distribution functions of the 
two ions are indistinguishable (dash dotted lines in Fig.[S]) 
and the center of the peak is located in between the two 
peaks of Ru^+ and Ru^+ in the charge constrained ET 
complex, at 2.14 A. This degeneracy is due to the elec- 
tron delocalization error of the BLYP exchange correla- 
tion functional. 



3 1 1 1 1 1 1 1 1 r 




Simulation time [ps] 



FIG. 6: Fluctuations of the vertical energy gap AE for the 
Ru'^'''-Ru^^ electron self-exchange reaction in aqueous solu- 
tion. The energy gap was calculated according to Eq. [14] for 
configurations generated with CDFT-MD. See section Hill for 
details. 

The energy gap Eq. [TJ] computed for an ensemble of 
configurations taken from the CDFT-MD trajectory, is 
shown in Fig. [S] (see section IIIII for computational de- 
tails). The average is {AE)a ~ 1-53 eV and the mean 
square fluctuation is (JAE'^)^''^ = 0.41 eV. The error of 
the average due to the finite length of the trajectory is 
estimated to be ~ 0.2 eV. The probability distribution of 
the energy gap fluctuations and the corresponding dia- 
batic free energy profile Eq. [T5] are shown in Fig. \7\ Due 
to the linear free energy relation Eq. [T7] two segments of 
the free energy curve are obtained from the CDFT-MD 
equilibrium simulation, one for the equilibrium region, 
and one for high free energies at the equilibrium region of 
the product state (see also section InC] and Refs. [38Ll47j). 
The two segments fit well to a parabola with a correlation 
coefficient of 0.99983. This shows that the Ru2+-Ru3+ 
electron-self exchange is well described in the linear re- 
sponse approximation, which is an essential assumption 
in Marcus theory of electron transfer. 

Since for electron self-exchange the reorganization free 
energy equals the average energy gap (Eq. [TB)) . we ob- 
tain X^{AE)a = 1.53 eV from CDFT-MD. This value 
includes the reorganization of the two Ru-hexahydrates 
and of 51 water molecules solvating the electron transfer 
complex. Reorganization of higher solvation shells and 
of bulk solvent is missing. In previous work we have es- 
timated a correction term for reorganization free energy 
of the missing solvent by carrying out a series of classi- 
cal molecular dynamics simulations for different system 



sizes and extrapolating the reorganization free energy ob- 
tained to the limit of infinite dilution ^43i]. The correction 
term for the 63 water molecule system amounts to 0.09 
eV. Our estimate for reorganization free energy of the 
infinitely diluted system is then 1.53 -I- 0.09 = 1.62 eV. 
This estimate is within the range of values obtained in 
previous classical molecular dynamics simulations of the 
same reaction using polarizable water models[43] , 1 .60 eV 
for SWM4-NDP water, 1.71 eV for P0L3 water and 1.87 
eV for AMOEBA water (values taken from Ref. [11], cor- 
rected for finite size but not for nuclear quantum effects). 
In a continuum studv[49j a value of 1.95 eV is reported 
that fits well the experimental rate constant. [5^ How- 
ever, this value was calculated for a Ru-Ru distance of 
6.5 A and should be corrected to 1.75 eV if the same 
distance as in the present CDFT-MD simulations is as- 
sumed (5.5 A). Unfortunately a direct experimental esti- 
mate for reorganization free energy is not available since 
the work term and the electronic coupling matrix element 
are unknown [50|. Yet, under a number of assumptions, 
a value of about 2.0 eV was reported in Ref. [50| . 



T 1 1 1 1 1 1 1 1 — I r 




AE [eV] 

FIG. 7: Probability distribution of the energy gap Eq.[l4](A) 
and diabatic free energy curves Eq.[T5](B) for the Ru^^-Ru^^ 
electron self-exchange in aqueous solution. Data points are 
collected in bins of width 0.11 eV. The data points in (A) are 
fit to a Gaussian and reflected about the origin to generate the 
distribution of the symmetrical product state (dashed lines). 
In (B) data points within 1.5 standard deviations from the 
center of the two distributions are fit to a parabola. 
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V. CONCLUSIONS AND OUTLOOK 

In this work we have presented an implementation of 
charge constrained density functional molecular dynam- 
ics in the plane- wave code CPMD. Several technical is- 
sues of CDFT were investigated such as the dependence 
of the results on the shape of the weight function used to 
define the charge constraint. Although for small donor 
acceptor-distances the energy depends strongly on the 
weight function, for medium to large distances the de- 
pendence is rather small. Thus, it is appealing to use 
this method for the study of long-range electron transfer 
problems. 

We have demonstrated that it is feasible to calculate 
electron transfer properties of condensed phase systems 
within the framework of CDFT-MD. As yet this is not 
possible with conventional density functional molecular 
dynamics simulation because for most donor-acceptor 
systems uncorrected GGA or hybrid density functionals 
give charge delocalized states for any nuclear configura- 
tion. Although the computational cost of CDFT-MD is 
higher than for standard Born-Oppenheimer dynamics, 
by a factor of about 2 — 3, CDFT-MD proved to be a vi- 
able method for sampling charge- localized diabatic states 
of condensed phase systems. 

The reorganization free energy obtained for aqueous 
Ru^+-Ru'^"'" electron self-exchange, 1.62 eV, is smaller 
than estimates reported previously by other authors, 
l.gSeVjiil and 2.0eV[53. This can be partly explained 
by the short Ru-Ru distance of 5.5 A in our present sim- 
ulation. Increase of this distance to 6.5 A0, [5l| will 
lead to an increase in reorganization free energy by about 
0.2 eV, thus bringing our value closer to previous esti- 
mates. However, present CDFT-MD simulations predict 
that the difference in Ru-0 bond lengths in Ru^"'"(H20)6 
and Ru'^"'"(Il20)6 is significantly reduced in the electron 
transfer complex relative to the isolated aqua ions. Thus, 
estimation of the inner-sphere contribution from ab-initio 
calculations [i^ or bond length differences [sol [sij of the 
separated aqua-ions could lead to an overestimation in 
reorganization free energy, because the deformation of 
bond lengths in the ET complex are neglected. 

In the present work we have focused on the free energy 
contribution to the electron transfer rate constant. Fol- 
lowing Van Voorhis and coworkers [2^ we will implement 
an approximate calculation of the electronic transition 
matrix element. This will enable us to predict absolute 
rates for electron or hole transfer in extended systems 
such as solvated donor-acceptor complexes and solids. 



APPENDIX A: WEIGHT FUNCTION CUTOFF 

The weight function defined in Eq. [5] is always finite 
even at points in space where all pi are zero. Yet, for nu- 
merical calculations we introduce for each atom species k 
a radial cutoff Rck to avoid small numbers in the denom- 
inator. This also makes the real space charge integration 



more efficient. The effect of the cutoff can be formally 
described by a Heaviside step function 9. Thus, in Eq. [5] 
we make the following replacement, 

P,{\y - R,|) ^ p,(|r - R,|)0(i?cfc - |r - R.|) (Al) 

where is the position of atom i. Rck is chosen such 
that the total reference density of species k is smaller 
than lO^^e, unless stated otherwise. 



APPENDIX B: CONSTRAINT FORCES 

The additional force on atom i due to the constraint 
term on the right hand side of Eq. [2] is given by 

Fc. = -l^/p(r)^^.r. (Bl) 

For the weight function defined in Eq. [12] the derivative 
in the integrand of Eq. IBll is given by 

g^^- /'('--^'D G.(r-R.) (B2) 
where 

{w(r - R,) - 1 ieD 

w(r - R,) + 1 leA (B3) 

w{r - R,) i(^DyjA 

depending on whether atom i is in the group of donor 
(D) or acceptor (A) atoms, or in neither of them (e.g. 
solvent atoms). The derivative of the density is given by 

P.:(ir-R.i) = 

^ gp.(|r-R.|) r-R, 

9|r-R,| |r-R,|- ^ ' 

The radial partial derivative of pi is calculated numer- 
ically using splines. However, due to the radial cutoff 
of the density (Eq. lAip the derivative Eq. IB4I has to be 
replaced as follows, 

pK|r-R.|) - p',{\v-nmRck-\v-m 

+ p,(|r-R,|)(5(i?^,fe-|r-R,|), (B5) 

where 5 is the Dirac delta function. Thus, the constraint 
force, Eq. IBll is composed of two terms, the force due to 
w within i?c/c, F^"'''^'', and a surface term, F™''^, 

F,, = F'.^'^^ + F^,^'-f (B6) 

where 
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Tiinsidc 



nsurf 



-y|p(r) ^"!j^^^^ e{R,k - |r - R.I) dv (B7) 



-Vpi{Rck)Rck 



p{Rck,^,^)G{Rck,^,^) 



Y,Pt{Rck,'&,f) 

sin d cos ip 

xR 

ck I sin d sin | sin i9 d-ddcp (B8) 
cos-i? 

The derivative in the integrand of Eq. IB7I is given by 
Eqs. IB2lB4l and in Eq. IB8I we have changed to spherical 
coordinates. The surface force Eq. IBSl is integrated over 
a thin shell on a cartesian grid. Rck times the radial unit 
vector is just the position vector of a point on the surface 
(x,y,z). Thus the surface element can be expressed in 
cartesian coordinates as 



sin?? d-ddip — sgn(z) 



y dxdz — X dydz 



Rck{Rck 



Z2) 



(B9) 



APPENDIX C: PREDICTION OF THE 
LAGRANGE MULTIPLIER V IN CDFT-MD 



for prediction of V{t) from the history of V. This should 
provide a better initial guess for the search of the un- 
known Lagrange multiplier. We chose a Lagrange poly- 
nomial of order n for this purpose [52|: 



k 

3=k-n 



n 



ti tj 



(CI) 



where k is the index of the last timestep. In Eq. ICll in- 
formation from n timesteps preceding step k is used to 
extrapolate the value of V for timestep k+1. Naively, one 
would expect higher order polynomials to perform better. 
Yet we found for CDFT-MD simulation of the aqueous 
Ru2+-Ru'^+ complex (section FlV Ep an optimum extrapo- 
lation order of n = 2. Higher order polynomials led to an 
oscillatory behaviour of Eq. ICll and to poor initial guess 
values. We note that the optimal order for extrapolation 
of V depends very much on the system under consid- 
eration. Therefore it seems advisable to calculate the 
optimum interpolation order for each system in advance 
from a short test-run before carrying out simulations on 
a larger scale. 
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